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Abstract 

We give a means for measuring the equation of evolution of a complex scalar field 
that is known to obey an otherwise unspecified (2+l)-dimensional dissipative non- 
linear parabolic differential equation, given field moduli over three closely-spaced 
planes. The formalism is tested by recovering nonlinear interactions and the asso- 
ciated equation of motion from simulated data for a range of (2+l)-dimensional 
nonlinear systems, including those which exhibit spontaneous symmetry breaking. 
The technique is of broad applicability, being able to infer a wide class of partial dif- 
ferential equations, which govern systems ranging from nonlinear optics to quantum 
fluids. 
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The laws of physics profoundly simplify an otherwise bewildering array of data 
collected by our measuring instruments [1]. Many of these laws take the form 
of partial differential equations governing the spatial and temporal evolution 
of field quantities associated with radiation or matter. These equations involve 
quantities that are often either single- component or multi-component complex 
fields. One may enquire whether it is possible to uniquely determine the asso- 
ciated equation of evolution from measurements made on such fields. In this 
context, there has been much research on such "identification" [2] of partial 
differential equations, both linear and nonlinear, given measurements of the 
field itself (see e.g., [3,4,5]). However, this earlier work restricts consideration 
to intrinsically real field quantities which are sufficiently slowly varying that 
both their amplitude and phase may be directly measured, providing data that 
can subsequently be used to determine the equation governing the evolution 
of the field. 
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Such existing methods are not applicable to complex fields whose phase is 
not known, whether these be intrinsically complex fields or rapidly-oscillating 
real fields described by a complex analytic signal [6]. In the case of complex 
quantum-mechanical fields, the phase is not directly measurable - however, 
following a line of investigation which dates at least as far back as Pauli [7], 
methods do exist for inferring wavefunction phase from non-interferometric 
measurements of probability density [8,9,10,11]. Regarding the case of classi- 
cal fields such as scalar electromagnetic disturbances, limitations of present 
detector technology imply that the phase of the disturbance is not directly 
measurable, at optical and higher frequencies [12]; following Gabor, such real 
fields are conveniently described using their complex analytic signal, which 
is so called on account of the fact that extension to complex time yields a 
function of a complex variable that is analytic in the lower half of the complex 
plane [6]. 

In the present paper, we demonstrate the feasibility of determining the par- 
tial differential equations that govern the evolution of complex fields, given 
modulus data alone. This comprises a means for "measuring" certain partial 
differential equations, which govern complex fields whose modulus is known 
but whose phase is unknown. 

Consider a complex (2+l)-dimensional scalar field governed by an equation 
belonging to the following class of dissipative, nonlinear parabolic equations 



Here, \I/ = ^f(x, y, z), (x, y) are Cartesian spatial coordinates, z is an evolution 
parameter such as time or propagation distance, Vj_ = (d/dx,d/dy) is the 
gradient operator in the xy plane, a is a real number, and both /(|^|) and 
gd^l) are real functions of a real variable. These last two functions respectively 
represent any non-dissipative and dissipative non-linearities which contribute 
to the equation of evolution. Special cases of the above class of equations 
are used to model a wide variety of both classical and quantum systems, 
such as monoenergetic electron beams [9], beamlike monochromatic scalar 
electromagnetic waves [13], intense scalar electromagnetic fields in nonlinear 
media [14], uncharged superfluids [15] and (2+l)-dimensional Bose-Einstein 
condensates [16]. An important subclass of equation (1) are those for which 
/(|^|) has its minimum value when ^ 0; in this case the field ^(x, y,z) 
displays spontaneous symmetry breaking [15]. 

When studying the applicability of a particular form of Eq. (1) to a given 
physical system, one often encounters the following logic, (i) A particular 
functional form of the equation is postulated or derived; (ii) this equation is 
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then used to predict the result of a given experiment or series of experiments; 
(iii) goodness-of-flt between the prediction and the data, often within the con- 
text of a parameterized model of the particular experiment being undertaken, 
is used as an indicator of whether the given equation is satisfactory - if the 
equation is unsatisfactory, one modifies step (i) and repeats the cycle until 
satisfactory results are obtained [17]. 

Notwithstanding the many successes of such a methodology, one may enquire 
if the process can be systematized. Specifically, can one determine (or "mea- 
sure") the equation of motion governing a complex field when the only 
measurable quantity is |\&|? In a novel approach, given in this paper, we give 
an affirmative answer to this question. 

A key assumption in our approach is that the desired equation belongs to 
the class of equations (1). We start by obtaining a "hydrodynamic" formu- 
lation of this equation via the Madelung transformation [18], ^(x, y,z) = 

\J I(x, y, z) exp[i&(x, y, z)\ , where I(x,y,z) = \^(x, y, z)\ 2 is the probability 
density (or intensity, as appropriate) of the field and y,z) is its phase. 
Substitute this into Eq. (1), and then separate real and imaginary parts, to 
obtain: 



f^ + V ± -(/V ± <t>) + Ig(I) = 0, (2) 

a Jz~ ~ /(/) + 1 V±$|2 " H{1, V±/) = °' (3) 

where the "diffraction term" H(I,'V±I) = /~2V^vT is a measurable func- 
tion of probability density or intensity (we will henceforth speak of intensity, 
rather than probability density, for convenience). The above pair of coupled 
equations is equivalent to the partial differential equation (1) from which it 
was derived. Note that, for matter- wave fields, the diffraction term disappears 
in the classical limit; for radiation wave-fields, this same term disappears in 
the geometric-optics limit. For either case, the classical limit of Eq. (3) yields 
the Hamilton- Jacobi equation governing the evolution of the field. 

The continuity equation (2) is independent of the functional form of the non- 
dissipative nonlinearity /(/). Introducing the scaled phase $ via $ = 2$/a, 
and restricting ourselves for the moment to the g(I) = (non-dissipative) 
case, this equation becomes: 

- + V ± ■ (IV ± $) = 0. (4) 



We now turn to the inverse problem of determining the phase of the field 
from intensity measurements alone. To this end we note the observation of 
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Teague, made in the more limited context of complex fields that obey the 
special case of Eq. (1) in which both / and g are equal to zero, that the 
associated continuity equation (4) may be used as a basis for deterministic 
phase retrieval [19]. Specifically if one measures both / and dl/dz over some 
plane z = zq, with / being strictly positive over some simply-connected region 
Q in that plane, and zero outside Q, then a simple non-linear generalization 
[10] of a result due to Gureyev and Nugent [20] shows that Eq. (4) can be 
uniquely solved for the scaled phase <&(x, y,z — z ), up to an unknown additive 
constant. We stress that this reconstruction requires knowledge of neither a 
nor /(/), relying rather on the hypothesis (which may be tested using the 
self-consistency argument given near the end of this paper) that \1/ obeys an 
otherwise unknown equation belonging to the class of equations (1). 

In retrieving the scaled phase <3>, one is faced with the determination of the 
appropriate boundary conditions for <3>. Here we briefly discuss the role of 
boundary conditions in the retrieval of $. To be specific we consider a phase 
problem in optics. As has already been stated, knowledge of both / and dl/dz 
over a plane z = zq, with I being strictly positive over some simply-connected 
region in that plane, and zero outside, implies that the continuity equation (4) 
(termed transport of intensity equation [19]) can be uniquely solved for the 
scaled phase <&(x, y,z — z ) up to an unknown additive constant [10,20]. When 
solving the transport of intensity equation for the phase, boundary conditions 
may be treated explicitly (see e.g., [21]), or implicitly, as in reference [20]. 
Regarding the explicit treatment of boundary conditions, let us give a non- 
linear generalization of an argument due to Roddier [22] (see also [21]), which 
shows how Neumann boundary conditions to the non-linear dissipative trans- 
port of intensity equation (2) may be determined from modulus data alone. 
In this example, z corresponds to a propagation distance, allowing us to place 
a thin black screen in the plane z = zo, into which is cut a simply-connected 
convex aperture with smooth boundary dQ. We assume the experiment to 
be so arranged that the intensity of the field, in the plane of the aperture, 
is strictly positive within the aperture, and zero outside. Since the aperture 
is both smooth and convex, points on the boundary of the aperture may be 
specified by the single- valued continuous function r = R(9), where r and 9 
are radial and angular plane polar coordinates with respect to a given origin 
which lies within dfl. The intensity distribution, in the plane of the aperture, 
can then be written as I(r,9, zq)H[R(9) — r], where H is the Heaviside step 
function, and I(r, 9, z ) is everywhere smooth and continuous. Substituting 
into Eq. (2), and writing in explicit functional dependencies for the sake of 
clarity, one obtains: 




H[R(9) - r] {V ± • [I(r, 9, ^ )V ± $(r, 9, z )\ + I(r, 9, z )g[I(r, 9, z )}} . (5) 
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Here, S is the Dirac delta and d/dn denotes differentiation along the outward- 
pointing normal to the aperture edge. Evaluate the above expression at the 
edge of the aperture; assuming the phase to be single-valued and continuous 
within and on the aperture, and the intensity to be differentiable within the 
aperture, the last term on the right side of the above expression are negligible 
compared to the first term. This leaves: 



\ ) z=zn 



(6) 



This equation may be radially integrated, to yield the Neumann boundary 
conditions which are required for the unique solution of Eq. (2) for the phase, 
up to an arbitrary additive constant. From a physical point of view, we note 
that the singular behaviour of the intensity derivative at the edge of the sharp 
aperture is due to local energy flow transverse to the aperture edge, as the 
field propagates from plane to plane. Note, also, that the above method for 
determining Neumann boundary conditions requires knowledge of neither the 
non-dissipative non-linearity / nor the dissipative non-linearity g. 

Evidently, knowledge of $ allows us to infer |V^$| 2 in Eq. (3). Further, <3> 
can (up to an unknown additive constant which, in general, differs from plane 
to plane) be determined over two closely-spaced planes z = z and z = 3z , 
so that one may infer d<&/dz (on the z = 2z plane) up to an unknown 
constant, denoted by R(z). Thus the estimated ^-derivative of $ is d<p/dz = 
d&/dz — R(2z ), so that Eq. (3) becomes: 



= 0. 



(7) 



If we assume that the intensity is non-vanishing in Q, but vanishes on its 

boundary, continuity of the intensity implies that at least two points 

and (x2,y2) i n ^ can have the same intensity. For such pairs, I\ = I2, where 



I x 



31 



Vj) 



Ij,j = 1,2. Evaluate Eq. (7) at each point, and subtract the 



resulting two equations, to give: 



712 - a 2 5u 



0, 



(8) 



where the known quantities 712 and 5 12 are: 



712 = 2 [H{h, V ± h) - H{I 2 , V ± I 2 )\ 

2 



0(0! -0 2 ) 1 

512 = dz + 2 



V ± $ 2 



(9) 
(10) 
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Note that the unknown constant R has been eliminated. Equation (8) gives 
ot = \J 112/^12, where a positive root is chosen because, as can be easily shown 
by applying the divergence theorem to Eq. (2) and then invoking conservation 
of integrated intensity, a > 0. Having inferred a, one can now infer the phase 
y,z = 2z ) of the wave-field via y,z = 2z ) = ^a<&(x,y,z = 2z ). 
Consequently, and bearing in mind that \I/ is completely specified by its mea- 
sured modulus and inferred phase, we are able to reconstruct even though 
a and /(/) (in the g = case of the equation of motion (1)) are unknown. 

Having solved the inverse problem of reconstructing the wave-field \I/ from 
measurements of its intensity, we turn to a second inverse problem which 
forms the core subject of the present note, namely a means for inferring non- 
linear parabolic equations of motion that govern a complex scalar field, given 
modulus data alone. For the case where g(I) = 0, the only unknowns are /(/) 
and an auxiliary constant R. Substitute the expression for a into Eq. (7) and 
rearrange, to obtain 



With the exception of R, each term on the right hand side is a known function 
of x and y. Therefore, the functional form of / can only be measured up to 
an unknown constant. In many practical cases, one knows a priori the value 
of /(0), in which case / can be fully determined. Note that the constant com- 
ponent of / may always be removed, via the transformation ^ — > ^ exp(inz), 
for suitable k. 

Thus far we have considered the g(I) — case of Eq. (1), which amounts to 
ignoring dissipation. We now show one means by which dissipation can be 
included in the analysis. Suppose that one prepares a z-directed plane-wave 
state, for which $ is independent of x and y, over the plane z = 2z . For 
such a state, which may have intensity modulations over the plane z = 2z , 
V ± • (A7_l$) = 0. Thus Eq. (2) becomes \adl/dz + Ig(I) = 0, therefore 
g(I)/a = —\d\vLl jdz. If one repeats this measurement for a series of values 
of J, one can determine g(I)/a, as a function of /. Having determined this 
function, multiply both sides of Eq. (2) by 2/ a, and make use of the fact that 
$ = 2$/a, to arrive at dI/dz+V±-(IV±$)+2Ig(I)/a = 0. Since the function 
g(I)/a is now known, and both / and dl/dz can be measured for a non plane- 
wave, the only unknown in this equation is <3>. This can be determined using 
the method cited earlier in the text, replacing dl/dz with dl/dz + 2Ig(I)/a. 
Having so determined both a and f(I) can be determined, using the same 
method as previously outlined. 

We now give a numerical example, to provide a simple demonstration of the 
application of our methodology using simulated data. As an initial condition 



6 



for numerical modelling using Eq. (1), we chose a zero-phase modulated Gaus- 
sian wave-packet centered on (x,y) = (x ,y ): 



Here, A and W are real parameters respectively denoting the peak squared 



and 5, n are real parameters. For the simulations presented here, we chose 
A = 10, W = 8, 5 = 0.01 and n = 20; the position of the peak is located at 
r = 0.5. 

In modelling the forward evolution of ^ from z = to z > 0, we must specify 
all terms in Eq. (1). We choose a = 1 for all simulations, which can always be 
arranged if one rescales z appropriately. The grid size is 2025 x 2025 pixels, 
with a z-step Az = 10 -7 , run for 5 x 10 3 time-steps. The squared modulus 
of ^ (intensity) is measured every 100 time-steps, so that at the end of the 
simulations, we obtain 50 measurements of the intensity profile. The intensity 
is measured over the central 1025 x 1025 sub-image of the larger simulation 
grid, corresponding to a spatial domain [x] x [y] of [0, 1] x [0, 1]. The spatial step 
of the simulation is therefore Ah = 1/1024. The evolution of ^(x : y, z — 0) to 
z > is obtained by finite difference discretisation of with d^/dz being 
approximated using fourth-order Runge-Kutta differentiation. 

The phase retrieval step, which involves solving Eq. (4) for $(x,y,z') given 
both / and dl/dz at a given value of z — z\ is amenable to the full multi-grid 
method [23]. Once $ and the derivative of $ are obtained from three slices 
of the measured intensity profile, we calculate a. For an intensity profile with 
maximum I m , we infer a at 0.2/ m , 0.4J m , 0.6J m and 0.8/ m . For each chosen 
intensity, we go through every point and check if the chosen intensity 

lies between the intensity at and (i + 1, j) or and + 1). If so, 
we store the point that has the closest intensity to the chosen intensity. We 
use the previously-described method to calculate a from each pair of points, 
and thereby construct a histogram of a. The value of a is determined from 
the peak of the histogram, with the error a a given by its full width at half 
maximum. 

Once a has been determined, we calculate $ from the definition $ = 2$/«. 
The result is then compared with the actual phase of the wave field \P taken 
from the simulations. As already stated, we can only measure phase up to 



y,z = 0) = ^/A(l + M x )(l + M y )e~^ r -^ 1 ) . 



(12) 



modulus and width of the packet, r = Jxl + yfc, 




M x = 5e 2l > w ) cos [2irn(x - x )] , 
My = <5e~2l-w^J s in [2nn(y - y )\ , 



(13) 
(14) 
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Table 1 

Measured values of a and /(/), for the first series of simulations, which have g(I) = 
0. The actual value of a is unity, with /(/) as described in the text. 07 denotes the 
RMS error in inferring /(/), after the constant offset has been adjusted (see text). 



a 


f(I) 




1.01 ±0.02 


-199.40 ± 100.26/ 


0.26% 


1.00 ±0.01 


-52.85 - 10.03/ 2 


0.28% 


1.00 ±0.02 


-297.17 + 401.19/ - 40. 12/ 2 


0.30% 


1.00 ±0.01 


-105.57 ± 20.30/ - 10.07/ 2 ± LOO/ 3 


0.66% 


1.01 ±0.04 


-118.43 + 100.01 sin(vr/) 


0.89% 



an arbitrary additive constant. For comparison, we shift the phase so that its 
average value is zero in a given plane, which corresponds to a particular global 
phase choice. The accuracy of the numerical method for phase recovery is 
dependent on the resolution of the numerical grid. It is found that the relative 
error in the phase, namely the normalised RMS error between the retrieved 
and the input phase, is of the order of 5% on average. It is expected that, 
since the multigrid method utilises smoothing and coarse grid correction, the 
major source of error arises from regions with large phase gradients. Further, 
since the phase is undefined when the intensity vanishes, phase retrieval is 
expected to be less accurate away from the interior region where the intensity 
approaches zero. 

Obtaining a and $ allows us to determine /(/) using Eq. (11). In the first series 
of simulations, which ignores dissipation, we used four different polynomials 
/(/) = 7 / + k/ 2 + C/ 3 with (7,k,C) = (100,0,0), (0,-10,0), (400,-40,0) and 
(20,-10,1), and a sinusoidal function /(/) = 100sin(7r/). The recovered a 
and /(/) are shown in Table 1. This demonstrates accurate recovery of both 
a and /(/) for nonlinear systems, including those with symmetry-breaking 
potentials. Measurement for the input /(/) = 100sin(7r/) (a = 1) is shown 
in Fig. 1. a can be read off from the peak of the histogram in (a). Accurate 
measurement of a leads to precise measurement of /(/) up to a constant as 
shown in (b). This constant arises as result of $ only being able to be retrieved 
up to an arbitrary constant, and this constant has been set to zero on both the 
first and the second planes. Since the constant part of /(/) can be transformed 
into the phase of ^ as pointed out earlier, it may be used to appropriately 
adjust the constant part of the recovered phase on the second plane. 

The second series of simulations incorporate dissipation, which was assumed 
to be of the form g(I) = DI@, where D and (3 are real numbers. The sim- 
ulations were performed for various values of D and f3 with a = 1 and 
/(/) = 10 sin(7r/). In this second series of simulations we show how dissipation 
can be measured, in addition to the other constants and functions which ap- 
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Table 2 

Input dissipation function, gi(I), compared to recovered dissipation g r (I)- {g r (I)} 
is the recovered dissipation, averaged over 49 measurements, while is the RMS 
error in this recovered dissipation. 

9i{I) 9r(I) <?g (9r(I)) °(g) 

10/ 9.43/ L04 2.11% 9.59/ 103 1.30% 
2/ 2 2.22/ 1 - 95 1.19% 2.00/ 2 - 00 0.94% 

pear in Eq. (1). Previously we pointed out that the nonlinear dissipation term 
g(I) can be measured using a plane wave formalism. For many systems it may 
be possible to construct a plane wave state for \1>; for example, if \1> describes 
monoenergetic electron beams or an electromagnetic wave field. However, for 
the systems where \1/ is a complex matter field describing an uncharged super- 
fluid or a Bose-Einstein condensate, it may be difficult to construct a plane 
wave for In this latter case the nonlinear dissipation g(I) may be mea- 
sured as follows. If the fluctuations of the system are small, the second term 
Vj_ • (JV_l$) in Eq. (2) is expected to be much smaller than the first and last 
terms. In this case the nonlinear dissipation term can be approximated by 

»C"-5S- (15) 



For a system with large fluctuations, we obtain the nonlinear dissipation by av- 
eraging the dissipation over many measurements. For N measurements, Eq. (2) 
can be written as 



N 

£ 

k=i 



adl „ 

2 8~ Z +V " 



(JV ± $) + Ig(I) 



0, 



(16) 



J k 



where k denotes the kih measurement in the ^-direction. The first and third 
terms in the square brackets in Eq. (16) are expected to increase with N, 
whereas the second (fluctuation) term is expected to approach zero in the 
limit of large N (long term evolution), for a chaotic system. Therefore the 
nonlinear dissipation can be estimated using the average value 

a N 

/ g m\ ~ -_V 

k=l 



ldl 

Idz 



(17) 



Note also that when measuring the nonlinear dissipation, the parameter a is 
not yet known. The nonlinear dissipation is in fact measured in terms of a. 
This is not an issue since g(I)/a is the only quantity needed to measure other 
functions in the evolution equation. For convenience we use the notation g(I) 
since we have set a = 1. The results for two typical simulations are shown in 
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Table 2. The fourth column shows a more accurate determination of g(I) when 
averaged over 49 measurements, which is an indication that the dissipations 
can be measured accurately for a system with large fluctuations (if a large 
number of measurements are obtained). Figure 2 shows the measured dissipa- 
tion functions for g^I) = 2I 2 and gi(I) = 10/. The dissipation function g(I) 
is accurately recovered, with an accuracy that can be improved by averaging 
over many measurements. 

The methods for measuring the evolution equation presented in this letter can 
be tested on empirical data; in this context we point out its constraints. To 
precisely infer the equation of motion from modulus information, we require 
high-resolution data. With the chosen grid resolution, our method does not 
cope well with Poisson noise much larger than one part in ten thousand. 
This currently restricts its application to high resolution and low-noise data. 
Notwithstanding this, there are currently many systems for which our method 
is applicable, such as the optics of intense electromagnetic beams in nonlinear 
media, and uncharged superfluid systems. 

We conclude by showing how the methods of this paper may be extended to 
multi-component (2+l)-dimensional complex fields, denoted by {ty n (x, y, z)}, 
which comprise a set of iV complex scalar wavefunctions \P n = ^ n (x, y,z),n — 
1, • • • , N. This multi-component wavefunction might obey a system of coupled 
nonlinear dissipative parabolic equations such as: 

d \ 
ia n — + V 2 ± + f n + ig n ]* n = 0, (18) 



where a n are real numbers, while f n ■ ■ • , In) and g n ■ ■ • , ijv) are real 
functions of N real variables, and /„ = |^/ n | 2 - The "hydrodynamic" formula- 
tion of Eq. (18) is: 



^ + V ± .(J B VJU + ^ = 0, (19) 

oz v ' a n 

Hn + fn-^ (^ + ^ V ^| 2 ) =0 > ( 20 ) 

where H n = J~ 1 / 2 V5_V^n is the generalized diffraction term. Equation (19), 
which is uncoupled and linear at the level of the unknown wavefunction phases 
$ n = |a: n $ n , is independent of /„. When g n = 0, this equation possesses a 
unique solution for the phases <3> n (each up to an additive constant), given I n 
and dl n /dz as data, provided all such phases are continuous. To subsequently 
solve Eq. (20) for a n and /„, our technique is to find pairs of points with 
the same /„. For arbitrary multi-component fields it is not known how such 
pairs of points can be found (although such pairs of points exist since f n 
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vanishes on the boundary and is non- vanishing in the interior). However, for 
two- component fields (N = 2) in two spatial dimensions, we can always find 
such pairs of points. That is we can construct a closed trajectory, T n , every 
point of which has the same intensity I n (n = 1 say). As we traverse a path 
in Ti, I 2 traverses the corresponding path in T 2 . However, since Ti is a closed 
trajectory, T 2 is necessarily a closed trajectory. If these trajectories possess 
any intersection points, then I\ = I 2 at such points. It is then straightforward 
to follow the methods developed in this paper to infer the equation of motion 
of the two-component field. A similar argument shows that in three spatial 
dimensions, it is always possible to infer the equation of motion of a three- 
component field. 

In order to measure equations of motion using the methods outlined in this 
paper, we require that the modulus of each field component be strictly positive 
over a simply-connected region Q. This implies that these wavefunctions can- 
not possess topological defects [20,15]. However, once the equation of motion 
has been inferred using a wavefunction free of topological defects, wavefunc- 
tion reconstruction in the presence of defects may be performed using the 
method outlined by Tan et al. [11]. 

Finally, we emphasize that our method employs highly redundant systems 
of equations, which yield multiple determinations for the desired equation of 
motion. If the postulated class of equations is insufficiently large, this will be 
manifest as a lack of internal consistency in the reconstructed equation, which 
thereby constitutes a testable hypothesis rather than an assumption. 
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(a) (b) 

Fig. 1. Measurements of (a) a and (b) /(/) at the end of a simulation (based on the 
input a = 1, /(/) = sin(7r7) and g{I) = 0). The histogram of a shows a sharp peak, 
and this peak is considered to be the measured value (denoted as a p ). Comparison 
with the input function /(/) = sin(-7r/) shows that we have accurately measured 
/(I) up to a constant of —118.43. This constant part is due to the fact that we 
only measure the phase of up to an arbitrary constant. Therefore this constant 
value may be used to correct our measurement of the phase of the wavefunction 
(see text). 




Fig. 2. Measuring the nonlinear dissipation for the functions g(I) = 10/ and 
g(I) = 2I 2 (a = 1 and /(/) = sin(7r/)). Plots with more oscillations are obtained 
from a single measurement (from data at the end of the simulation), whereas the 
less oscillatory lines are taken from averaging over 49 measurements. This illustrates 
that the dissipation can be measured more accurately by averaging over larger sam- 
ples. Averaging over many measurements may be a viable method for measuring 
nonlinear dissipations for highly fluctuating systems. 
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